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Abstract 



We present a closed analytical formula for the scattering intensity from 
charged hard sphere fluids with any arbitrary number of components. Our 
result is an extension to ionic systems of Vrij's analogous expression for un- 
charged hard sphere mixtures. Use is made of Baxter's factor correlation 
functions within the Mean Spherical Approximation (MSA). The polydis- 
perse case of an infinite number of species with a continuous distribution of 
hard sphere diameters and charges is also considered. As an important by- 
product of our investigation, we present some properties of a particular kind 
of matrices (sum of the identity matrix with a dyadic matrix) appearing in 
the solution of the MSA integral equations for both uncharged and charged 
hard sphere mixtures. This analysis provides a general framework to deal 
with a wide class of MSA solutions having dyadic structure and allows an 
easy extension of our formula for the scattering intensity to different poten- 
tial models. Finally, the relevance of our results for the interpretation of small 
angle neutron scattering experimental data is briefly discussed. 



I. INTRODUCTION 



Most of the naturally occurring or industrially important fluids are actually mixtures 
of more than two components. In particular, liquid mixtures of macroparticles (e.g., mi- 
cellar or colloidal suspensions) exhibit size, shape and possibly charge distributions of their 
components. Macroparticles, unlike atoms, are intrinsecally polydisperse. 

As a consequence, any rigorous statistical mechanical model of real mixtures should be 
able to include a large finite number p of species, which eventually will become infinite in 
the presence of polydispersity. In the latter case the properties of the mixtures will depend 
upon idealized continuous distributions of size, charge, etc.. 

Unfortunately, the presence of a large number of components poses challenging problems 
to the two major tools of the liquid state theory, namely computer simulation and integral 
equation (IE) methods. In fact, in the first case a prohibitively large number of particles 
may be required for a satisfactory Monte Carlo or molecular dynamics simulation, with 
consequently long and expensive computational times. In the second case, when the lEs 
must be solved numerically (as almost always occurs), the algorithms work with arrays 
formed by about 10^ x elements, so that larger and larger amounts of memory and again 
of computational time are required as p increases; moreover, non-convergence problems in 
the numerical solution are encountered more frequently, especially in the presence of strong 
correlations and high density. Finally, polydisperse systems with infinitely-many components 
cannot be investigated in a fully adequate way if the lEs have to be solved numerically. 
In such discretization is unavoidable and one is forced to replace the continuous 

distribution of species with an appropriate p-component mixture, as shown by D'Aguanno 
et al. 0, 0], 0], p who solved integral equations with p < 10 for Yukawa (= screened 
Coulomb) plus hard core interaction potentials. However, apart from the aforesaid work, 
most of the published IE studies on fluid mixtures concern binary or, more rarely, ternary 
systems [!,§,§. 

In this paper we are interested in the static structural properties of multicomponent 
fluids. Experimentally, these can be efficiently monitored using X-ray, neutron and light 
scattering. Fitting experimental structure factors or scattering intensity of systems with a 
large number of components using IE theories is actually feasible only if the relevant equa- 
tions do admit analytical solutions. Such a fortunate opportunity occurs only for peculiar 
interaction potentials and, to the best of our knowledge, this is limited to a single class of 
approximate equations. 

The simplest case refers to the model mixture with hard sphere (HS) interactions. Using 
Baxter's analytical solution of the approximate Percus-Yevick (PY) integral equation for 
HS mixtures 0, Vrij []TU| derived a closed expression for the scattering intensity R{k) which 
holds true for any number p of components. It is remarkable that Vrij's formula depends on 
p only through the averages of some quantities. Hence, only the number of terms involved 
in these averages increases as p increases. This feature allows to extend the application of 
this equation for R{k) to polydisperse HS mixtures [0. The results are in a satisfactory 



agreement with Monte Carlo simulation data [IT2 



Independently and along a somewhat different route, Blum and Stell [1^ presented the 
PY scattering function for polydisperse fluids of hard or permeable spheres in terms of 
Fourier transforms Hij{k) of the total correlation functions (some misprints are corrected in 
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Appendix A of ref. [jT4[). In addition, the averages involved in Blum and Stell's expressions 



were evaluated analytically by Griffith et al. [T^ in the case of polydisperse mixtures with 
a continuous Schultz distribution of diameters. 

Now, apart from the HS case, we are unaware of a similar analytic determination of 
scattering functions for any other interaction potentials ||I6[, |T^. The main purpose of the 
present paper is then to extend Vrij's work to multicomponent mixtures of charged particles. 

The simplest model for charged systems is the so-called primitive model, consisting of 
an electroneutral mixture of charged hard spheres embedded in a dielectric continuum and 
interacting through unscreened Coulomb potentials. For these pair potentials, the unre- 
stricted general solution, within the Mean Spherical Approximation (MSA), is available ||18|| , 
T9| , [pO| , ||2l[| : simple analytical expressions for the Baxter factor correlation functions qij{r) 
IP are given and only a single non-linear equation in one parameter, F, has to be solved. 

It is well known that, for dilute solutions of highly charged particles, the MSA may 
lead to unphysical negative values for some of the radial distribution functions gij{r) close 
to the contact distance or in a neighborhood of the ffist minimum (in fact, the MSA is 
asymptotically correct for r — > oo, but it is incorrect at short distances, in the region just 
outside the core). To overcome this drawback and extend the validity of the MSA to arbitrary 
low densities, some authors have proposed an ad hoc rescaling method, which replaces large 
electrostatic repulsions by appropriate HS repulsions |Q, . Such a "rescaled mean 



spherical approximation" (RMSA) preserves the analytical form of the MSA solution and 
therefore allows an iterative fit of low density experimental data. 

Nevertheless, in the regime of weakly charged solutions at sufficiently high concentrations, 
the MSA is expected to be a reasonably accurate approximation (in fact, the Coulomb part 
of the potentials can be considered as a perturbation with respect to the HS one, whereas 
in the above-mentioned opposite regime the electrostatic effects predominate over the HS 
repulsions). 

The primitive model can be utilized for a large class of ionic fluid mixtures (electrolyte 
solutions, molten salts, solutions of macromolecules and polymers, micellar and colloidal 
suspensions, microemulsions, etc.) and with a careful analysis of the regime of validity for 
each case the MSA can be safely employed. For instance, Abramo et al. ||2^ applied the 



MSA analytic solution for charged HS to the evaluation of partial structure factors and 
x-ray diffiaction patterns for the whole family of molten alkali halides (binary systems). 
Caccamo and Malescio ||2^ compared MSA and HNC (hypernetted-chain) results for struc- 
tural and thermodynamic properties of polyelectrolytes (in particular, micellar solutions 
with an added electrolyte). We also note that Senatore and Blum already performed 
MSA calculations of the average structure factor S'^'^lk) for mixtures of charged HS with 
either size polydispersity or charge polydispersity. However, none of these authors presented 
closed- form analytical expressions for the MSA structure functions as in Vrij's work [|lOl; this 
is the task we have accomplished here. 

The paper is organized as follows. In the next Section we briefly review the basic for- 
malism of scattering and integral equation theory for multicomponent fluids. In Section HI 
we discuss some useful properties of a peculiar class of matrices (related to dyadic matrices) 
encountered in the MSA solution for both neutral and charged HS mixtures (for uncharged 
HS the MSA coincides with the FY approximation). Section IV then provides the analytical 
MSA equations required to derive the new formula for the scattering intensity (and for the 
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average structure factor). This is described later on in Section V. The polydisperse hmit 
is briefly discussed (Section VI) and, finally, in Section VII we summarize our results and 
comment on perspectives for future investigation. 



II. SCATTERING AND INTEGRAL EQUATION THEORY 



A. Scattering functions 

The coherent scattering intensity I{k) for a p-component fluid mixture with spherically 
symmetric interparticle interactions can be written in terms of the partial structure factors 
S{j (^) as [ p8|| 

Rik) =I{k)/V=J2 imf^ F^ (k) F, (k) S,, (k) = (k) S (k) f {k) (1) 

Here, k is the magnitude of the scattering vector, V is the volume of the system, the 
number density of species i. Further, Fj(fc) = F^Bi {k) denotes the scattering amplitude (or 
form factor) of species z, with being the scattering amplitude at zero angle and Bi (k) the 
(angular averaged) intraparticle interference factor. The components of the column vector 
f (k) and of its transpose, the row vector {k), are the form factors Fi [k) weighted by 
means of the corresponding densities. More precisely, we define 

m = p]"F,{k) (2) 

In order to get more compact formulas, often our notation slightly departs from Vrij's one 
|T^. Specifically, here this author uses fi to denote our 

Finally, S {k) is a symmetric matrix whose elements are the Ashcroft-Langreth partial 
structure factors Sij (k), defined by 

Sij (k) = 6ij + Hij (k) = 5ij + {piPjf^ hij (k) , (3) 

or, more concisely, 

S(A:) = I + H(A;), (4) 

with I being the unit matrix of order p {6ij = Kronecker delta). Here, hij (k) is the three- 
dimensional Fourier transform of the total correlation function, hij (r) = Qij {r) — l, and Qij (r) 
is the radial distribution function between two particles of species i and j at a distance r. 

In the short- wavelength limit (A; ^ oo) , when the particles scatter the incident radiation 
independently, one has S (A; — oo) = I and therefore 

R{k^oo) = j: p,Ff (k) = (k) f (k) = |f ik)\' (5) 

i=l 
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where |f denotes the magnitude of f (A;). The normahzed scattering intensity is also 
called the effective structure factor or the measured average structure factor 

c^nu^- _ F(A:)S(fc)f(A:) _^ , F (fc) H (fc) f (fc) 

^ ^ = /(A:-oo)- ^(^yp ^(^yp 



B. Direct correlation functions 

A second possible representation of the scattering functions can be obtained in terms of 
the three-dimensional Fourier transforms Cij {k) of the direct correlation functions Cij (r) . 
These functions are defined by the Ornstein-Zernike (OZ) integral equations of the liquid 
state theory, which, for systems with spherically symmetric interactions, read as 

P r 

hij (r) = dj (r) +^ Pm j dr' dm {r') /i^^- (|r - r'|) (7) 

m=l 

Note that these equations can be solved only if they are coupled with a second relationship 
between cij (r) and hij (r) . Such a "closure" consists of the exact formula 

(r) = exp [-(3(Pij (r) + 7^^- (r) + Bi^ (r)] - 1, (8) 

plus an approximation to the "bridge" function Bij (r) , which is a complicated functional 
of hij (r) and higher order correlation functions [^] ( (pij (r) is the interparticle potential, 
/3 = 1/ (ksT) , ks is Boltzmann constant and T the absolute temperature; jij (r) = hij (r) — 
Cij{r) ). 

By Fourier transforming the OZ convolution equations, these can be written in fc-space 

as 

[I + H (k)] [I - C (k)] = I, (9) 
with Cij (k) = {piPj)^^"^ Cij (k) . Since S (A;) = I + H (A;) , we also get 

S (k) = [I - C {k)]-\ (10) 

or, equivalent ly, 

^^Ak)- \l_C(k)\ ' ^^^^ 

where |I — C {k)\ is the p x p determinant of the matrix 1 — C {k) and |I — C (A;)P* is the 
cofactor of its (j, i)th element. 

By substituting eq. ( p!0D into eq. (|^, we then obtain the scattering intensity in terms of 
the Cij (fc), i.e., 

i?(A;)=f^(A:)[I-C(A:)]-^f(A;), (12) 
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while the average structure factor is always given by S^{k) = R{k)/ |f {k)f . 

For theoretical investigations all these formulas may be more useful than those in terms of 
the hij (k) . In fact, in the cases in which the OZ equations have been analytically solved (for 
some particular interparticle potentials and with appropriate "closures") the (r) have, 
in general, rather simple expressions, whereas the hij (r) do not admit a simple analj^ic 
representation and are usually evaluated by numerical inverse Fourier transform of h^j [k) 
(in fc-space, the hij (k) are again much more involved than the Cij (k) ). 



C. Baxter factor correlation functions 



Vrij |T0[ proposed a third representation of the scattering functions, in terms of the so- 
called Baxter factor correlation functions qij (r) . Baxter showed, for hard sphere fluids, 
that the OZ equations can be transformed in an equivalent, but often easier to solve, form, 
by introducing a Wiener-Hopf factorization of the matrix I — C{k) . Later on Hiroike ||3^ 
extended Baxter's work to disordered fluids with any kind of spherically symmetric potentials 
and obtained generalized Baxter equations without using the Wiener-Hopf factorization. 

Noting that I— C(fc) is a symmetric matrix and an even function of k, Baxter |^ suggested 
the following factorization 



I-C(A;) =Q^(-/c)Q(A;) 
where the elements of Q (k) are of the form 



(13) 



Q (A;) = I - Q (A;) , with Q (k) 



dr e'^'Q (r) 



(14) 



and 



Qij{r) = 27T {pipjf^ qij{r) 



(15) 



( Q"^ is the transpose of Q ). In general, in the analytically solvable cases, the factor 
correlation functions %(r) and Qij{k) have a even simpler mathematical form than Cij{r) 
and Cij{k), respectively. 

Substituting eq. (|13D into eq. (|10D yields 



S(A;) = Q-^ (k) \Q-' i-k) 



(16) 



and, since 



(k) 



Q{k) 



Qik) 



(17) 



then 



S^Ak)=J:Q^r^ik)QJ^i-k) 



D{k) 



Y.\^{k) Qi-k) 



mj 



(18) 
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with 



D{k) = Q{k) Qi-k) = |I-C(A;)| 



(19) 



For the scattering intensity, plugging eq. (|T6D into eq. (|l]) and assuming that f (k) 
f (— k) yields the simple expression 



R{k) = (k) s (-k) = —y^L^ (k) L (-k) 



Dfk) 



(20) 



where we have defined 



sfk) = 



Q-(k)l"f(k).iM 
Q(k) 



(21) 



The hypothesis f (k) = f (— k) , which is equivalent to Fi (k) = Fi {—k) , is really correct for 

the physical systems we are concerned with in this paper. In fact, for homogeneous spheres 
with radius o"j/2, it results that 

Fi (k) oc Vi^, ikai/2) , (22) 

where Vi = {tt/6) af is the volume of a particle of species i and $i (x) is an even function 
which may be expressed in terms of the first-order spherical Bessel function ji (x) as 

$1 (x) = 3ji (x) /x = 3 (sin x — x cos x) / x^ (23) 

Eq.s (1^) and (pOD give the partial structure factors and the scattering intensity in terms 
of the Baxter functions Qij{k). 



III. THE ANALYTICALLY TRACTABLE CASE OF DYADIC MATRICES 

In view of producing theoretical expressions for the scattering functions, we have pre- 
sented three possible routes, based upon different correlation functions, i.e., Hij {k) , Cij {k) 
and Qij{k), respectively. 

In Baxter's route, followed in this paper, the crucial point for the analytic evaluation of 
the scattering functions is the possibility of getting a simple and closed expression for the 
inverse (k) of the matrix Q (k) . We want to emphasize that such a task can easily be 
accomplished if Qij{k) has a particularly convenient structure of the following general form 

n 

Qvik) = + a^t^hf^ = 5,, + % (k) (24) 

Clearly, Wij (k) is simply related to the unidimensional Fourier transform qij (k) of qij (r) , 
i.e., Wij (k) = -2tt {pipjf^"^ % (k) . 

It is now convenient to recall that a second-rank tensor Tij = aibj formed by the outer 
product of two vectors is sometimes represented by a symbol called a dyad, ab, and a linear 



7 



combination of dyads A^a^^^b^^) is called a dyadic [Q. In addition, as we propose here, 
a dyadic formed by a sum of n dyads might be called an n-dyadic. 

Using this terminology, we could say that eq. (|2^) is equivalent to require that qij (k) is 
a generic n-dyadic matrix (of order p) . 

In the case of the PY solution for neutral HS mixtures (including the polydisperse limit), 
Qij {k) is a 2-dyadic matrix [^, |T^. In the MSA solution for charged HS |19|, as we shall 



better see in Section V, qij (k) becomes a 3-dyadic. 

Both Vrij and Blum with coworkers [l^ investigated some properties of these 
matrices, independently and from complementary points of view. On the one hand, Vrij 
emphasized the role of the rank of Q (k) and developed an original method to derive Q^^ (k) . 
On the other hand, Blum and coworkers gave more directly closed expressions for {k) for 
both 2- and 3-dyadics (on passing, we note that none of these authors uses our terminology; 
for instance, Blum and H0ye [|19| refer to these matrices as Jacobi matrices, whose inverse 
can be found for any size of the matrix). 

In order to join and generalize the two above-mentioned approaches, in this Section we 
shall present some general properties of n-dyadic matrices, which hold true for any n (the 
relevant proofs are outlined in Appendix A). 

Let us start from the fact that for any matrix of the form 



6ii + W, 



(25) 



where Wij is a generic p x p matrix (i.e., not necessarily a dyadic), the determinant |M| can 
be expressed by the following expansion in terms of principal minors 



IMI 



1 + W 



1 ^ 
3' ^ 



P 1 P 

1 + + 



i=l 



w,, W,, 
% W,, Wjk 

Wk^ W,, Wkk 



+ 



33 



(26) 



For a generic Wij matrix, this expansion clearly stops at the p-th term. However, it is 
easy to show that, if Wij is an n-dyadic (with n < p), then any its minor having order m > n 
vanishes! Using this along with some manipulations, we find our first basic result on dyadic 
matrices: any p x p matrix of the form 



M = I + a^'^^b^'^) 
11=1 



(27) 



always has rank n, irrespective of p. Consequently, the expansion ( |^^ of its determinant 
breaks off after the nth term and \M.\,which is of order p, turns out to be equal to a determi- 
nant Dm of order n < p,as expressed synthetically by 



IMI 



Dm = det ((5^, + a(^') ■ b^'^) 



(28) 



where det(...)^ means the determinant of a matrix of order n, while the dot denotes the 
usual scalar product of vectors 
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(29) 



It is to be noted that the use of expansion ( ]26|) with 2-dyadics is originally due to Vrij 
T0| ; the extension to n > 2 and eq. (|28|) are new. 

Let us now consider the cofactors |M|"'* which will be needed in the calculation of the 
inverse matrix Q"^ (k) and, consequently, of the scattering intensity, as we shall discuss in 
Section V. In the same Appendix A it is also shown that for a matrix of the type p5| ) with 
a generic Wij the cofactors can be written as follows 



IM 



Ml 6ij - U, 



where the following expansion for Uij holds true 



E 

k=l 



Wi 



kj 



Wkk 



+ 



E 

k,l=l 



+ 



(30) 



(31) 



Again for n-dyadic matrices this expansion breaks off at the level n < p, as it is straight- 
forward to check. Along the same lines as before, we find that Uij can be written as follows 



(32) 



where the D'f^ are n x n determinants which can be obtained from the determinant D 



M 



upon replacing the fi-th. row with a row formed by bj 



(1) 



in) 



From M-j^ = |M|^V |M| , it is now immediate to get the inverse matrix M 



E« 



(m) 



D 



M 



with 



D 



ip) 



det ( [1 - 5^^] + a(°) ■ b(^) + 5^^ 



(33) 



(34) 



This general formula constitutes our second result on dyadic matrices, yielding a new 
simple and elegant recipe which holds true for any n value: the inverse of any p x p matrix 
M = I + X;5^=i a(^)b(^) has elements given by eq. as linear functions of simple deter- 

minants of order n.The determinant is obtained from D]\j, defined in eq. (|^^, by 

using a Cramer-like rule, i.e., by replacing the fi-th row with a row formed by bf \ b^^ (of 
course, all rows can be interchanged with the corresponding columns, without altering the 
value of the determinant). 

Due to the symmetric role played by the a's and the 6's in eq. (0), it is clear that 
one could derive an alternative expression to (|32| ) with the role of the the a's and the 6's 
interchanged. 

Note that our formulas for |M|''* and are fully new with respect to Vrij's work []10 . 
Eq. (33) agrees with the expressions of Mj~^ given by Blum and coworkers [13|, |19] for 
n = 2 and n = 3; at the same time, our generalization offers a systematic way to calculate 
M^-*- for any n. 



We shall be using all these results in Section V. 



9 



IV. CHARGED HARD SPHERES AND MSA SOLUTION 



The primitive model for ionic fluids consists of a p-component electroneutral mixture of 
charged hard spheres embedded in a continuum of dielectric constant e (which may represent 
a possible solvent). The species i has diameter at, number density pi and electric charge ZiC 
(e is the elementary charge). The interparticle potential is the hard sphere Coulombic one, 
i.e., 

+00 for r < a-ij = | (cTj + aj) 

cP,, (r) = <; (35) 



e ZiZj/{er) for r > a 



and the electroneutrality condition reads 

PmZm = (36) 



m 



In r-space the Baxter form of the OZ equations which can be derived from eq. (|13]) is 



' rcij (|r|) = -gjj(r) + 2tt pm IZ^ 9™ (t) qj^j (r + t), r> \j = \{ai- aj 



(37) 



. rhij (|r|) = -q/Jr) + 27r Em pm I^^ dt qim (t) (r - t) h^j {\r - t\) , r> X 



where the prime denotes differentiation with respect to r. Note that, if Ajj < 0, then r 
may assume negative values as well, although the correlation functions Cij{r) and hij (r) are 
determined, of course, only for positive distances. 

The Mean Spherical Approximation, used to solve these equations, consists of adding to 
the exact hard core condition 

hij (r) = — 1 for r < CTj^, (38) 



the approximate relationship 



Cij{r) = -p(pij{r) 



for 



r > a, 



(39) 



which is asymptotically correct for r 
with the PY approximation: (r) = 



00. For uncharged HS potentials, the MSA coincides 



{exp 



(r)] . The MSA, 



eq. 



may also be regarded as the following approximation to the bridge functions 



B,, (r) = In [1 + 7,,(r) - /50,,(r)] - [7,,(r) - /?0,,(r)] 



for 



r > a 



(40) 



The solution of the MSA integral equations for the most general case of ions with arbi- 
trary charges and diameters was found by Blum |18|, |19], |2^, [^| 



It is worth noting that 

Blum employed a form of the Baxter equations which is slightly different from eq.s (P7[). 
In this paper it proves convenient to follow the original form proposed by Baxter [^j and 
Hiroike |3^, eq.s (|37D , and used by Vrij ||10[. The differences from Blum's formulas of ref. 
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T9| are, however, small (a detailed comparison is given in Appendix B). Our modified form 
of the MSA solution reads 



fo, 



r < Xij 



|(r - aij)^ g/' + (r - aij) q/j - AiZj, Xij < r < a. 



r > aij 



(41) 



The coefficients of the factor correlation functions are given by 

1 / 36 \ r2 



-AiAj 



(42) 



(43) 



with 



' r 1 + r^i 



6 = ^Ep™< (^ = 2,3) 



(44) 



(45) 



A = 1 - 6 



?B = Pe^/e (Bjerrum length) 



(46) 
(47) 

(48) 



2 e 1 + 



(49) 



Finally, the value of the parameter F can be determined by solving numerically the 
following consistency equation 



(50) 



which is equivalent to the condition 



(51) 
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In the limit of point ions (all am ^ 0), 2r tends to the Debye inverse shielding length 
kd of the classical Debye-Hiickel theory for electrolyte solutions, while for finite size ions 
(2r)~^ is always larger than : 



1/2 

(52) 



There is a striking similarity of the MSA to the Debye-Hiickel theory, with 2r taking the 
place of kd (2r is the correct screening parameter for finite size ions). 
It is to be noted that the degree of the algebraic equation (|50D for T increases with the 
number p of components in the mixture. Moreover, among all its solutions, only one is 
physically acceptable: the one which is positive and tendes to kd from below in the infinite 
dilution limit. 

V. APPLICATION OF DYADIC PROPERTIES TO CHARGED HARD SPHERES 

A. The inverse matrix (A;) 

The starting point is the evaluation of the unidimensional Fourier transform, qij (k), of 
the MSA solution for qij (r). The resulting Qij (k) can be written as 

Q,, (k) = 6,, + % (k) = 6,, + {p,p,f' W,, (k) , (53) 

with 

Wi, {k) = -27cq,, {k) = e"""^/^ [aj {k) + a,^, (k) + A^j {k)] 

(54) 

= Zi, {k) + e^^'A,7, {k) 
where i - when it is not a subscript - is the imaginary unit and 

Xm = kaj2 (55) 

47r 1 vr 1 

= ^^'""^^ ~ = 6^^' A*^^ ^^'^ ^^^^ 

f5,{k) = f3o,j{k)+Pi,j{k) (57) 
^ = f f ^ r^'l*^" ^^^^ ^^^^ 

7, (k) = 70, , (fc) + 71, , (A;) + 72, (fc) (60) 
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70, j {k) 

71, 3 (^) 

72, j (k) 



P. A 



2i 



aj (k) + -^/?o, j (k) 



(61) 
(62) 
(63) 



Here, the function $i (x) = 3ji(x)/a; is the same as in eq. (|23D and we can also write 
$0 (x) in terms of Bessel functions: $o (x) = 3jo{x) = 3sinx/x. 

Note that aj {k) and /5o, j {k) are even functions, unhke [3j {—k) = (3* {k) and 7^ {—k) = 
7* {k) (the asterisk denotes complex conjugation). Hence, we have Wij {—k) = W*j {k) and 
4, {-k) = % (k) . 

Moreover, Zij (k) is the pure hard sphere term considered by Vrij, while e'^^^Ai'^j (k) is 



the electrostatic contribution (a^ and Pj correspond to Mj and A'^- of ref. [|T^, respectively). 

In the limit {all charges Zm} — > 0, Ai,Pz and 7j(A;) vanish simultaneously. In the long 
wavelength limit (fc 0), Zij{k) tends to a finite value, whereas 'jjik) diverges as k~^. Such 
a singularity of qij {k) implies a k~'^ divergence of Cij (/c); however, all this does not infiuence 
the structure factors, which remain finite at the origin. 

The crucial point is now that, according to eq. (|5^) , the MSA expression of W^j {k) 
for charged HS mixtures is a 3-dyadic. As we discussed in Sec. II, this peculiar form of 
Qij (k) allows to write down a simple and compact expression for its p x p determinant 



Q{k) 



I-Q (k) 

gii(A:) • 



Using the result 
■ QiM {k) 



^ we get immediately, for any p. 



DQ{k) 



Qmi (k) ■■■ Qmm (k) 
where we have defined 

{fg) = T.Prr^'^'^'^fm{k)9m{k) 



(«) iP) (7) 

{aa) l + {a(3) (^7) 
{Aa) (AP) 1 + (A7) 



(64) 



(65) 



Notice that our shorthand notation is somewhat different from Vrij's one. Of course, (a), 
{13), etc. are complex functions of k. However, as in eq. (p^), we shall mostly omit the 
argument, unless when necessary. The connection with the a's and 6's is 



(66) 





- ^1/2 iX. 
Hi ^ 1 


if 


1/2 

= p/ aj 




1/2 iX 




= Pj 


(3) _ 
i 


„l/2 iX, 4 . 

Pi e /ii , 




1/2 

= P/ 7i 



Eq. ( p^ clearly stresses the fact that, even if the order of the determinant Q (k) may 
become quite large with increasing the number p of components in the mixture, nevertheless 
Dqik) remains of order 3 {contraction or reduction effect). 
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From 



Q (fc) = Dqik) and eq. (|T9|), it follows that 

D{k) = |I - C (k)| = DQ{k)DQ{-k) = DQ{k)D*{k) 



(67) 



Once again, D{k) = |I — C (k)| is a determinant of order p, whereas DQ{k)DQ{—k) is a 
product of two determinants of order three. Since |A| |B| = |AB| (if A and B are any two 
square matrices of the same order), DQ{k)DQ{—k) could be rewritten as a single determinant 
of order 3. We can thus conclude that, under the peculiar assumptions of our model, also 
the matrix I — C (k) has always rank 3, no matter how large is the number p of components. 

Similar considerations hold true for the cofactor matrices. From eq.s (|30|) and (|3^) one 
can easily show, after some manipulations, that 



Qik)''= Qik) 5^, - ip^p.y^'' e'""' Df\k) + a:,Df\k) + A,Df\k) 



{hi = i,--,p) 

(68) 



where Dj^^ (/i = 1,2,3) is simply defined as the 3x3 determinant derived from Dq by 
replacing its mth row with the vector {aj, Pj, 7j) . 

Then, from Q"-^ (k) = Q{k) / Q{k) , it is now immediate to get the inverse matrix 

of Q(A;) 



Q-i (k) = 1-V{k) 



with 



Vij{k) = {pipjf'^Vij {k) = (piPjY 



a'{k) + ai(3'{k) + Aa'{k) 



(69) 



(70) 



a'{k) = 



Df\k) 
DQ{k) ■ 



3 '(k) = ^ 



(71) 



Here, the prime does not denote a derivative: we have used this notation only to emphasize 
the similarity of the form of Vij{k) to that of Wij{k) in Qij {k). 



B. Scattering intensity 

To get the scattering intensity 

R{k) = Y.s,{k) s,{-k), 



(72) 



where 



Sj{k) 



kLj {k) 



(73) 



let us first evaluate 
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jm 



(74) 



Note that in eq. ([73|) both the numerator and the denominator have been multiphed by k 
to "heal" the k~^ singularity of 7j(fc) at the origin and ensure that in numerical calculations 
R{k) has a regular behaviour as k —>■ 0. 



If the expression 



for 



- - , , J"* 

Q [k) is inserted in this equation, then, after some manip- 



ulations, Lj [k) can be cast in a very convenient form as a 4 x 4 determinant, i.e.. 



Lj (k) = pj 



1/2 



(F) 


1+ (a) 


(5, 




(7) 


(aF) 


(era) 


1 + {a(3) 


((T7) 


(AF) 


(Aa) 




1 + (A7) 



(75) 



However, to perform the summation over j, it is more convenient to re-expand this 
determinant along the first row 



L, (k) = p'/' [F,{kmk) + a,{k)T2{k) + f3j{kmk) + ij{k)T,{k)] 



(76) 



where T^{k) {p = 1,...,4) is the cofactor of the (l,/i)th element of the first row. Clearly, 
T,{k) = DQik). 

In this way we obtain 



R{k) 



with 



f +{Fa}Cl + {F7} Q 

+ +{a^}CiC*^ +{af3*}CiC; + {a^*} CiC; 

+ {PF}C2 +{/?«} C2Q +{Pr}C2CI 

[+{lF}C, +{^a}C,Ct +{^f3*}C,C; +{^Y}CsC; 

{fg} = J2 Pmfm (k) Qm (k) 

m 

C,{k)^T,^,{k)/T,{k) (/i = l,2,3) 



(77) 



(78) 
(79) 



Clearly, {-F^} is identical to f^f. 

This simple and elegant expression for R{k) arises in a quite general way from the dyadic 
structure of qij{k) (only some properties of a-,- were added, i.e., aj{—k) = aj{k) = a*{k) ). 
Of course, such a formula might be used as is. Nevertheless, from a practical point of view, 
it is more convenient to take advantage of the particular, explicit, expressions for aj^ j3j 
and 7j in the present model, in order to achieve a sort of "reduction" of the above-written 
result. Since simpler expressions are preferable for numerical calculations, we have tried to 
perform all possible simplifications, starting from the determinant form of Lj (k) and using 
some properties of determinants. 
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The first step is to take the determinant form of Lj (k), eq. (|75D, and subtract from the 
third column the second one multiplied by (8^2/^ — ik/2) . Since j is given by eq. (^9|), 
such on operation "eliminates" j from the third column, which becomes 

(/5o) -'-t + ¥k 

(80) 

1 + w 

The advantage of the resulting expression for Lj {k) is that, in place of the complex functions 
it involves only the f3o, j, which are real and even. It is worth noting that such a reduction 
is implicit in Vrij's work |I0 |. 



Let us now transform the electrostatic contributions. Our second step "eliminates" 72, j 
from the fourth column. Looking at eq. (p5|), it is evident that it is now necessary to subract 
a linear combination of the second column and the third one (in the new form with /5o, j), 
multiplied respectively by Pz A and 2iPz A /k. The fourth column becomes 

701, j = 70, j + 71, j 
(701) + P^2Pz 

(81) 

(0-701) - f P.A 
1 + (^7oi) 

As a further step, we have observed that the particular functional form of 71, j{k) allows 
a direct analytic evaluation of the necessary averages relevant to this term. By deriving from 
eq.s (H) , (^8l) , (^) and from the electroneutrality condition the following new identities 

{A} = Y^p^A^ = -^Pz (36 + LA) (82) 

rr, Vrl 



{a A} = J2 PmamAm = A, (83) 
together with eq. ( |^ ) for {A^} = J^mPmA"^ = ^s/vr , we find 

(71) = -t^Pz (36 + TA) (84) 

(a7i) = t'^PzA (85) 

2r 

(A71) = ^^ (86) 
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Putting all these results together and multiplying Lj{k) by k to "regularize" its behaviour 
at k = (as proposed in eq. (|73D we get 



kLj(k) = p 



1/2 



(F) 






^7oi, j 


1+ (a) 


(/?o)-36/A + zfc/2 


(7o) - 27riP^rA 


(aF) 


{pa) 


1 + {o%) 


(a7o) 


(AF) 


{Aa) 


(^/?o) 


A; [1 + (A70)] + ilY 



(87) 



where 



701, =70, j(fc)+7i, 



TAj- 27r 
77 T 



2r 

sin ■ + i ( cos X,- H — — sin X, 



(88) 



This new form for kLj{k) is really more convenient than that corresponding to eq. ([75D, 
since it requires a smaller number of numerical calculations and allows to write a simpler 
expression for the scattering intensity. In fact, expanding the aforesaid determinant along 
the first line and taking into account that Fj, aj and /3o, j are real and even functions, we 
get after some algebra the following final result for R{k) 



R{k) = Ri{k) + R2{k), 



(89) 



R^{k) = {f'} + CiCl + {l3l} C2CI 

+2Re \{Fa) Ci + {F/5o} C2 + {a/3o} CiC*] 



(90) 



R2{k) 



+2k Re [{7oii^} C3 + {701a} + {7oi/3o} C3C*] 



(91) 



where Re [...] denotes the real part of a complex value function and the used here are not 
the same as those in eq. ([77|): their definition is again given by eq. ( |7DD but the cofactors 
T^{k) now refer to the new determinant expression for kLj{k), i.e., to eq. (|87|). Note that 
T,ik) = kDQ{k). 

Observe that {Fa} , {F^q} and {a/^o} are real quantities, while the 70, j{k) are odd and 
purely immaginary functions, i.e., 70, j{—k) = —70, j{k) = 7q^ ■{k). 

In the limit of {all charges Zm 0}, it results that R2{k) and Ri{k) — > R^^{k), 
i.e., we recover exactly (apart from the notation) Vrij's formula for the neutral HS case. 

Finally, as concerns the partial structure factors, they can easily be evaluated by means 
of eq.s (|16D and (|69|)-([7T|), but we prefer not to give their explicit expressions here. 



VI. THE POLYDISPERSE LIMIT 



In the treatment of polydisperse systems it is customary to start from a mixture with 
a finite number of components and then to perform a limiting process, which results into a 
continuous distribution of the properties which characterize the species {identity variables) 



34 1 . In other words, in an IE approach one first solves the equations and then takes the 



infinite-species limit. This is a delicate point. 
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Still in the IE framework, we like to emphasize a second alternative route, i.e., the 
possibility of taking the infinite-species limit from the outset. According to this second 
viewpoint; one first sets up integral equations for polydisperse fluids and then tries to solve 
them. It is worth mentioning that, although a polydisperse version of the OZ equations 
for HS has already appeared in the literature we are unaware of any papers where 



polydisperse Baxter integral equations have been proposed. 

To this aim, let us first recall how the composition of polydisperse systems can be speci- 
fied. In a "fully polydisperse" mixture the discrete variable 1 < i < p that labels the species 
is replaced by a set of one, two or more continuous variables, which constitute an identity 
vector I = {Ii, I2, Is} (not to be confused with the previously used identity matrix). The 
identity variables are distributed according to a multivariable molar fraction density func- 
tion p(I) [|^, [^: the molar fraction Xi = Pi/p = Pi/J2mPm of the discrete case goes 



mto 



p(l)dl = p(Ii, I2, Is)dlidl2...dl„ (92) 



which is the probability of finding a particle with species labels in the identity range 
(I, I-|-dI). As a consequence, for the number density it results that 

^ pp{l)dl, (93) 

with p being the total number density. 

In the opposite case of "fully monodisperse" mixtures, one could write 

p(I)=Ex(r)5(I-I'), (94) 

{!'} 

where {!'} is the set of identity vectors corresponding to discrete species present in the finite 
amounts {x (I')} (here, 6 is the Dirac delta distribution). 

Of course, one may also consider "partly polydisperse" mixtures in which several species 
are monodisperse, while the remaining ones are polydisperse. In this case, the distribution 
function p(I) is the sum of a discrete part, PdiJ), like that in eq. (0) and a continuous 
one, Pc(J)- Clearly, to be a probability density function, p(I) must satisfy the normalization 
condition 

^x(I')+ /f/Ipc(I) = l (95) 

{!'} 

Bearing all this in mind, it is immediate to pass from the Baxter equations (0) for 
p-component mixtures to those for polydisperse fluids, using the following the replacement 

Y,Pm ■■■ / dpm ■■■ = P dl p(I) (96) 

if we re-interpret the species label m as a shorthand notation for the identity vector I, whose 
distribution function p(I) may include both a continuous and a discrete part, respectively 
for polydisperse and monodisperse components. 

Using this simple convention, the procedure of solving the polydisperse Baxter equations 
for uncharged or charges HS, within the MSA, is formally the same as in the monodisperse 
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case, apart from the substitution required by eq. (0). In other words, we can conclude 
that, at least for the two aforesaid systems, the MSA solution for the polydisperse case can 
be obtained from the monodisperse one via the simple prescription expressed by eq. (pG]). 
Note that all this fully agrees with the findings by Stell and coworkers [|l^], |^ for neutral 
HS in the PY approximation and with the MSA investigation performed by Senatore and 
Blum on mixtures of charged HS with either size polydispersity (and fixed charge value) 
or charge polydispersity (and fixed diameter). 

Before concluding this Section, it is to be noted that, while for HS a single identity 
variable - the diameter - is sufficient to characterize a species (i.e., I = {o"}), for charged 
HS two variables - the diameter and the charge - are necessary (i.e., I = {a, z}). To avoid 
the difficulty of working with a two- variable distribution p {a, z) , one might assume a strong 
correlation between a and z, reducing the number of independent variables to only one. It 
is customary to express such a relationship in the form z = z (cr) |2^, [Q], but perhaps the 
inverse function a = a {z) might be preferred in view of the integrations to be done over 
the remaining independent variable (recall that the range of a is [0, +oo[, while z varies in 
the interval ]— oo, +cx3[ ). Anyway, for macroions of colloidal suspensions it is reasonable 
to assume that their charges scale linearly with their surface area, i.e., z oc cr^, keeping the 
surface charge density constant ||^, 0|. 



VII. CONCLUDING REMARKS 



The first and main result of this paper is given by eq.s (|89D -(|9lD, which provide a new 
MSA formula for the scattering intensity from multicomponent mixtures of charged hard 
spheres. This fills an important gap present so far in the relevant literature. It can be applied 
to a wide variety of ionic fiuids and is expected to be particularly useful to fit experimental 
data. 

Our result extends Vrij's one for neutral HS and reduces to it in the limit of vanishing 
electrostatic contributions. Despite the presence of the electrostatic terms, this case is still 
analytically tractable, due to the dyadic form of qij{k). Indeed, the rank of the matrices 
involved in such calculations turns out to be independent of the number p of components 
in the mixture, always being equal to the number n (typically <C p) of dyads constituting 
the n-dyadics qij{k). Such a feature allows to take the polydisperse p oo limit with an 
arbitrary continuous distribution of both charges and sizes. 

Our findings are by no means academics. The importance of having an analytical ex- 
pression of the scattering intensity for experimental purposes has been well established |^ 



An application of the present results to neutron scattering experimental data will be the 
subject of a forthcoming paper. 

A word of caution on the range of applicability of our results is in order. 

The failure of the MSA closure at low densities and high charge values could be overcome 



by using a rescaling procedure (RMSA) [|^, p3|, However, this does not present a 

problem for many ionic fiuids of current interest, which are made up of weakly charged 
particles at high volume fractions. In such cases the MSA is expected to be a reasonably 
accurate approximation, yielding semiquantitatively good predictions for the static structure 
factors. 
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Regarding the non-linear equation for the single parameter F, which has to be solved in 
the MSA scheme, selecting the physical root from the manifold of solutions is not a terribly 
difficult task |3^. In fact, as Hiroike pointed out, this equation has only one real positive 
solution which is the physical one. Moreover, in their application of the MSA for charged 
HS to micellar solutions Senatore and Blum [^] already solved the equation for F even in 
the polydisperse limit without mentioning any particular problems. 

As far as the potential considered in this work and the possibility of extending our MSA 
scheme to other different interactions are concerned, it is evident that the primitive model 
represents the simplest physically significant choice to describe ionic fiuid mixtures. 

The primitive model includes both attractive and repulsive Coulomb potentials and has 
the advantage of treating large and small ions on the same footing. For colloidal suspensions, 
it definitely gives a better description than the essentially one-component models which 
consider only the repulsive (screened) Coulomb interactions between macroions, while the 
other species (counter-ions, any added electrolyte and solvent) are taken into account only 
upon determining the screening in the effective macroion-macroion potential ^8|, p2[, ||^. 



Of course, several alternative choices of potential model may be found in the literature 
on scattering from multicomponent fluids. For instance, very recently, Kline and Kaler 
fltted their experimental partial structure factors of a mixed colloidal system (sodium 
dodecyl sulfate micelles plus Ludox colloidal silica) by a model of hard and sticky hard 
spheres (in the terminology of these authors, "sticky" means that an attractive square well 
is added outside the hard core). 

Indubitably the next step after the primitive model should be the mixture of par- 
ticles interacting by HS plus Yukawa (= screened Coulomb) potentials, i.e., /30jj(r) = 
Kij exp{—nr)/r for r > Cij. Unfortunately, the general MSA solution for multicomponent 
Yukawa fluids |^T]], |2T|] is much more involved than the MSA one for unscreened charged 
HS, even if substantial simpliflcations occur for factorizable coupling parameters, i.e., when 
Kij = KiKj In particular, a satisfactory treatment of polydisperse Yukawa sys- 



tems is still lacking. In this context, we recall that Lowen et al. |^ proposed a mapping 
of the polydisperse Yukawa model onto the polydisperse HS reference system. It is also 
worth mentioning that, to avoid the difficulty of working with the MSA solution for Yukawa 
mixtures, Ruiz-Estrada et al. invented a method to map a Yukawa system onto an 
equivalent (but much easier to be treated) primitive model system. Very recently, an ana- 
lytical equation of state for the HS Yukawa polydisperse fiuid (in the equal-diameter case 
with polydispersity in the coupling parameters Ki) has been presented [^.We believe that 
the Yukawa model deserves further investigation, specially on the possibility of deriving a 
closed MSA formula for the scattering intensity in this case as well. We are planning to do 
it in the near future. 

Finally, before concluding this paper, we remark that our methodological recipe, based on 
the properties of the dyadic matrices presented in Section III, has far reaching consequences. 
It is indeed the property of "contraction" or "reduction" of these matrices (according to 



which 



evaluation o 



can be reduced to an equivalent determinant of very low order) that makes the 
' the structure factors for fiuids with a large or infinite number of components 
a tractable problem. Our findings have also led to a simple recipe for building the inverse 
matrix Q^^(A;) for the most general qij{k) of dyadic form. Our formula not only includes 
those previously found by Blum and coworkers for n = 2 and n = 3 [jl3[, |jT9|, but, more 
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importantly, it provides a systematic way to extend the results of the present paper to new 
more complex cases with n > 3. 
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APPENDIX A: PROOFS FOR SECTION III 



I) The proof of (|26D is patterned after a similar one used in the Laplace expansion of the 
characteristic polynomial of a p x p matrix. 
Consider a p x p matrix of the generic form: 



6ij + Wij 



6,, + M'/'Wij 



(Al) 



For the sake of simplicity, we shall work out the case p = 3 first and then show how the 
result can be generalized. Let us rewrite the determinant of the 3x3 matrix (Al) in such 
a way that each element is the sum of two terms (type A and type B respectively) 



|M| 



l + PiWn 
+ {p3PlY^'Ws^ 



+ {P1P2Y^'W^2 0+(piP3)'/Wi3 
I+P2W22 0+ (p2P3)'/'W^23 

+ {P3P2Y^^W32 I+P3W33 



(A2) 



Using elementary properties of determinants, the above expression can be written as a 
sum of a determinant containing only elements of type A (i.e., the identity matrix), three 
determinants containing two rows of type A and one of type B, three determinants containing 
one row of type A and two of type B and, finally, one determinant containing only rows of 
the type B. Hence there are 2^ determinants overall. Again exploiting elementary properties 
of determinants and symmetry considerations on the permutation of the indices, one easily 
establishes that (A2) can be written as 



|M| = l + ^PiWH + -E 

i=l ^- ij=l 



+ P1P2P3 



Wn W12 Wi3 
W21 W22 W23 
W31 W32 W33 



(A3) 



The above expression can be generalized to the general p-component case, by noting that 
all terms with s < p rows of type A and p — s oi type B can be written as 



pp- 



w. 



s+l,s+l 



w. 



s+l,p 



p,p 



(A4) 



Collecting terms of all orders yields then (|26|) . 

II) Let us now show that, if Wij is an n-dyadic (with n < p), then any minor of W 
having order m > n vanishes. 

We outline the proof only for n = 2, since its generalization to n > 2 is very easy. We 
have to show that any minor of order m > 3 is zero, but, since all minors of order m > 3 
can be expressed in terms of the third order ones, it is sufficient to demonstrate that 



If 1^, 



W,, W„ Wi, 

W,, Wjr Wjs 
Wk, Wkr Wks 



a^^h^P + af'^b^ and if we define 



\/{i,j,k,q,r,s) 



(A5) 
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Ai = af^/af^ and /x, = bf /hf (if af ^ 0, bf ^ 0), 



(A6) 



then we could write Wij = a'pb^p (1 + XifJ^j). Substituting this expression into the determi- 
nant of eq. (A5) and factoring out all common factors, we get 



ct^ (X j dj^ b^ bj bj^ 



1 + XiHq 1 + XiHr 1 + Xifis 
1 + Xjllg 1 + XjHr 1 + XjIXs 
1 + Afc/ig 1 + XkHr 1 + XkHs 



(A7) 



Then, by subtracting the third row from both the first and the second one, we see that the 
resulting determinant vanishes since two rows are proportional. 

Ill) As concerns |M| expressed in terms of + a^'^) ■ h^'^\ eq. (^81) can be verified very 
rapidly for n = 2, by using eq.s ( ]26|) and (A5) together with the following identity 



(1) (2) _ (2) (1) 



_ 5(2)5(1) 



(A8) 



which holds true for 2-dyadics. 

IV) Finally, in order to prove eq.s (|30|)-(|32D for the cofactors, one first establishes that 
the elements M^^ of the inverse matrix can be computed order by order in the p's, by using 
the definition 



(A9) 



fc=i 



This leads to the following result 



fc=i 



k,l=l 



(AlO) 



Then, starting from |M|j' = Mij"'"|M|, using the expansion (^) for |M| and collecting all 
terms of the same order, one finds, after some simple algebra, the required expansion. 



APPENDIX B: COMPARISON BETWEEN DIFFERENT NOTATIONS 

In place of Baxter's factorization expressed by eq. (|TB|), Blum and coworkers |T^, [|TH[, 
m, start from 



I-C(k) = Q(k)Q^(-k) 
where the elements of Q {k) are again of the form 

Q(fc)=I-Q(A;), 

but with 



(Bl) 



(B2) 
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Q (fc) = 27r / °° dr e^^^Q (r) and Qij(r) = (pipj)^/' qij(r) (B3) 

This choice leads to a shghtly different version of the equations for the factor correlation 
functions, namely to 

27rrcij (|r|) = -9ij(r) + Em Pm IZj ^jrn (t) Qi^ (r + t), r> Xji = \ {uj - a>) 

(B4) 

27rr/iy (|r|) = -q-j{r) + 27r Em Pm /a^^ o?^ - /i^m (k - t|) gmj (i) , r > A^^ 
The correspondence between our (or Baxter's) notation and Blum's one is very simple 

^Ba.ter^^^ = ""^(0 (B5) 

As a consequence, our g^j, g/' and are easily obtained from Blum's counterparts by 
performing index transposition and division by 27r. 

Finally, our and correspond, respectively, to A^'^P^ (or A~^P„) and r2A of Blum's 
papers. 
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